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Abstract —This flow in rivers is concerned with unsteady 
flow in open channel and it mathematically governed by the 
Saint Venant equation, using a four-point implicit finite 
difference scheme. For a one-dimensional applications, the 
relevant flow parameters are functions of time, and 
longitudinal positions. Considering the equations for the 
conservation of mass (continuity) and conservation of 
momentum. The mathematical method is empirical with the 
computer revolution, numerical methods are now effective 
to develop the hydraulic model. Conventional methods used 
to compute discharge, such as stage-discharge and stage- 
fall discharge relationships has been inadequate. Further 
research is recommended to include tributaries in rivers. 
Keyword—Finite difference scheme, Numerical methods, 
Momentum equations. 

I. INTRODUCTION 

For engineering purposes, most of the solutions of the 
unsteady flow equations are numerical and are determined 
by using digital computers. There have been a great number 
of different numerical techniques that could be used. Some 
can be discarded as being inaccurate or unstable or too time 
consuming. However, there is still no answer as to which is 
best. 

Further are two fundamental divisions of methods for the 
solution of hyperbolic partial differential equations, 
describing one-dimensional unsteady flow in open 
channels: The hypothesis of the one dimensional flow is 
very often problematic when rivers with flooded plains are 
in question, so that the question is often raised whether the 
mathematical model should be one-dimensional or two- 
dimensional, Doyle, W. H. Jr. (1983). In the case of two 
dimensional flow there are three independent variables: two 
spatial ones (x) and (y), and the time variable (t), assuming 
gradually varied flow conditions the flood propagation 
equations could be written by analogy to the tidal equations, 
Faye, R. E. (1984) . The system of three equations is 
analogous to the system derived by St. Venant for the flow 
in one spatial dimension. The main assumptions are the 
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uniform distristatic vertical pressure distribution, Brakke 
K. (2008) 

In engineering practice different methods are used for the 
two-dimensional modeling of flood-plains, the quasi two- 
dimensional or “cell” model, and the mathematical two- 
dimensional model, explicit or implicit. On the basic 
principles of cell model, the flood plain remains completely 
submerged only during the flood peak. Out of this period, 
the plain still forms a system of cells separated by dikes, 
irrigation channels, culverts, etc. all of them influencing the 
direction of the flow. Bakes, I. (2015). If it is possible to 
assume that the flood plains, as well as the river bed can be 
divided into a certain number of cells, then each cell 
communicates with the neighboring cells and the links 
between them correspond to an exchange of the flow, 
Xiaoyong Zhan (2003). The limits of the cells are defined 
by natural boundaries, such as are roads, dikes, natural bank 
levels, etc. 

The cell centers must be defined so that the direction of the 
flow is correctly allowed for at any time during the flood. 
The slope of the water surface is defined for each cell (i) by 
the corresponding horizontal levels (Zi), assumed at the 
center of the cell, Fread, D. L. (2008). The area ABCD can 
be treated as one cell, but in this case the model would not 
reproduce the slope between cross-section AB and CD. This 
means that the area ABCD must be divided in two cells: 
ABEF and EFCD, with representing water level at the 
centers of these cells, (T’) and (T”) respectively. For each 
cell (i) the continuity equation and discharge equation 
between this cell and the adjacent cells must be written. The 
continuity equation for a cell can be derived based on the 
following two fundamental hypotheses: The volume (Yfi of 
the water stored in cell (i) is directly related to the level 
(Zj), and the dischargew between two adjacent cells (i) and 
(k) is a function of the corresponding levels (Z t ) and (Z k ). 
The change in the volume of water stored in cell during 
time (At) can be defined. 

The method of characteristics deals with problems imposed 
by nature, but often with serious difficulties. However, this 
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method has a very important role in determining the 
boundary values in many fixed-point finite difference 
methods. 

The general aspects of the method of characteristics was 
explained for St. Venant's equations written in the form of 
qe (1) 

Momentum Equation 

^ = S 0 - S f (2) 

g dt g Ox B dx 1 

Combining the two basic equations with the two equations 
of the total derivatives of two dependent variables (V) and 
(A), 

r)ii r\ it 

(3) 


1 du + u 
dt g 


du 1 
dx B 


ixr du i . dv , 

dV = — . dt H- .dx 

dt dx 

i a dA , , dA j 

dA = — . dt -i - .dx 

dt dx 


(4) 


and by solving the system of these four equations, the 
resulting characteristic equations are obtained: 


dx I A 

— = V + a.- 

W*- - _ I D 


dv+ p- 

— A.B 


dA = 


(S 0 -Sf)-3±%-. 


(5) 


(6) 


The equations (5) - (6) are four ordinary partial differential 
equations with the unknown values being (x; t; v; A). The 
integrals are line integrals along the characteristics. No 
approximations have been made. 

A(x,t 0 ) — A{x),and (7) 

V(x,to)=V(x) (8) 

The problem has to be defined for t > and ( t 0 ) 
corresponds to the steady flow regime. 

We assumed that in points 1, 2, 3,..., see Fig. (IB), Eqs. (7 
and 8) that represent the initial conditions are known. 

The objective of this paper is to fill the ( x, t) - plane with 
the net of the characteristics, so that the dependent variables 
can be defined at the intersections. This means, that by 
starting from "known" points 1, 2, 3, ..., the values of two 
dependent and two independent variables in points 5, 6 and 
7, can be calculated by solving the system of equations (5) 
and (6). Following the above described procedure, the 
whole ( x, t) - plane is filled with characteristics, as shown 
in Fig. (1) so that the dependent variables are defined at the 
intersections. 

Stability is tested by asking if a particular part of the 
solution, perhaps an initially unimportant part, is likely to 
grow without limit until it destroys the calculation. 

The stability is defined by the following well known 


equation: 

At 1 

Ax “ |V 0 ±C 0 | 


(9) 

which is known as the Courant - Fewy - Friedrichs 
condition, or often simply the Courant condition, where 
(F 0 ) - is the initial value of the mean flow velocity; (C 0 ) - 


the initial value of the small wave celerity. The boundary 
conditions are specified at a boundary, the number of 
conditions must be exactly equal to the number of 
characteristics originating at that boundary. 

II. METHOD 

The four possibilities for external boundaries are shown in 
Fig. (2) 

In case a), the condition is specified at the upstream 
boundary and is used with the equation of the negative 
characteristic (R-M) to calculate the second dependent 
variable at point (M). 

In case b), one condition is given at the downstream 
boundary, and is used with the equation of positive 
characteristic along (F-M). 

In case c), the dependent variables must be specified along 
the boundary so that the values at (M) are known prior to 
the solution. 

In case d), no boundary values are specified and the 
unknowns must be interpolated between (F-M) (or R- M), 
to obtain the dependent variables on the boundary. 

In literature one can find many methods of characteristics, 
one of which will be presented as follows: 

The basic equations (5) and (6) will be written in the 
following form: 

- the pair of equations used for defining the positive 
characteristic: 

dx = (v + y/g.A/B). dt 

( 10 ) 

dF + + a + . dt = 0 
( 11 ) 

- the air of equations used for defining the negative 
characteristic: 

dx — v — 

( 12 ) 

dF~ + a~ . dt = 0 

(13) 

Where: 

F + + U + - 
— 2 

(14) 

U = —■ / ~7= 

(15) 


(16) 

If the derivatives in the equations (14-15) are replaced by 
finite differences, the basic equations can be written in the 
following form: 


±0-(S f -S o ) + fJg 
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%m ~ %u — a u- (/m ~ tu~) (17) 


F±-F u + a+. (t M - t u ) 

(18) 

Z D — f D . ( t M — tjf) 

(19) 

F m ~ F d + a D . (t M — t 0 ) = 0 

Where: 

(20) 

«U = (y + Jg-A/B) u 

(21) 

Pd = (V- Jg-A/B) D 

Other symbols are explained in Fig. (3) 

(22) 

By using the equations (17-22), two independent unknowns 

in point (M) are calculated first, and after that, 
(F^) and ( Fm) are defined. 

functions 

Two dependent variables in point (M) can be 
easily from: 

calculated 

4 = v f + u M 

(23) 

FM = - V f+U M 

(24) 

U M = U(Z) 

(25) 

THE FINITE DIFFERENCE METHODS ON A FIXED 

GRID 


In contrast to the characteristic methods, most hydraulic 
engineers dealing with river hydraulics prefer the fixed grid 
methods, in order to solve the equations of motion at a finite 
number of grid points in the ( x, t) - plane. There are two 
basic types of fixed grid finite difference schemes: the 
explicit and the implicit scheme. In explicit scheme, the 
equations are arranged so as to find the solution for one 
point at a time. In an implicit scheme, a group of advance 
points is solved by using simultaneous equations, which 
include the unknowns at all points in the group. 

III. RESULT 

Finite differences using a fixed rectangular net in the ( x, t) - 
plane but, it is very important to notice, that the theory of 
characteristics, even it is not used directly, comes into play 
in deciding the relative values of (A t) and (A x) which will 
ensure convergence. The criterion of stability is the same as 
for the method of characteristics, defined by Eq. 9, this 
being the most serious disadvantage of explicit methods. 
One of many existing variants of the explicit method of 
finite differences will be introduced, in order to explain the 
most important details of calculation. 


The derivatives are replaced by finite differences in the following way: 

dZ ~ Z(p,i+i)~Z(p,i+i) t dQ ~ Q(p,i+i')~Q(p,i+i) 
dx 2Ax ’ dx 2Ax 

dZ ~ Z(p,i+i) — ^(p,i+i) _ dQ ~ Q(p,i+i)~Q(p,i+i) 
dt At ’ dx Ax 

f(p,i+ 1) — 2 ' I f(p,i+i) + f(P.i+ 1)] 


Where: 

(Z)and ( Q ) the average values ; (/) — any function-, (p) — time row, (t) — distance row 
For the explanation of basic principles of this method, the St. Venant’s equations will be written in a form that corresponds more 
to the natural watercourses with a non-prismatic channel: 

- Equation of continuity: 


Where: 


dA dQ 

Tt + ^-~ qe 

—— . — ^ F — —-. —— T (1 — Ff) . — — — 

g.A dt g.A dx v ry dx 

(27) 

(28) 

= — 1^1 - — -S f +D e 

g.A \ dx\ g.A * e 

(29) 

(F r ) - the Froude number; 


V — W e cosa 

D e — --- ,q e ; (W e ) — velocity of the lateral flow, 

g. A 



(a) — angle of the lateral flow, 

\dA/dx\ z=const ~ the term involving the influence of the non 
The discretization according to Esq. (27), leads to: 

- Equation of continuity: 


prismatic channel. 


A (p+i,i) A (p,i± l] _j_ Q(p,i+1) Q(p,i+1) 
At 2 Ax 

- Momentum equation: 


"f Re(p,i+ 1) 


= 0 


(30) 
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Q(p,+i,i) Q(p,i± l) _|_ 2 -Q(p,i± l) 


At 


A(p,i± l) 


Q(p,i+i) Q(p,i-i) 


2 Ax 


_ Q \p,i± l) 
A 2 ,, 


\ A(PM1) A(pM1 A +g.A 

L 2 Ax \ 7=7n y 


(P.i± 1) 


] + [iM.p, ■ 1: 

[ 


Q (p,i±i). B (p,i±l) 


A 2 




j (p,i+1) z (p,i- 


2Ax 


.tl) j 


5 /( P ,i±l) +D1 (W± 1 )J 


Q (p,i±i) j[ 1 (p,i±i) 

^(P<i±i) 


(31) 


(p,i± l) L 2 - fl * J Z=Z 0 

Equations (30) and (31).represent a pair of algebraic equations, that can be used for the explicit determination of dependent 
variables Q + 1 ^ and Z( Pi+1 ^ in the unknown point (p+l,i), see Fig.(4). The above described scheme is valid for the so called 
"interior" points. For the "boundary" points (points 1. And D in Fig. 4), the upstream boundary condition and the negative 
characteristic are used for calculating the points in cross-section (1), while the downstream boundary condition and the positive 
characteristic are used for the calculation of the points in the downstream cross section (D). 

The equations used for calculating the points at the 
left boundary, in cross section X = Xi are: 

- the upstream boundary condition: 

Qu — Quit) — inflow hydrograph 

- the negative characteristic equation: 


Q(p+i,i) Q(p,i) 


■f 


At + [ 7 (p,D _ C (p,D 

z (p+i,i)~ z (p-i) , r I7 

ft LHp.U 


Q(p,2)-Q( p ,i) 


At 

Rp.i)~Ap. 2) 1 , r 

■[ Ax L =r + £?e (P-l)‘ L (P'l) 


Ax 

r ^ Z (p,2)~ Z tp. 

L (P, 1)J 


+ B 


(P, i) 


K. 


i) 


Ax 


-} = 9-\ v , «■[- 


Ap.DJ 


^/(p, i) 


+ D 


e(p,1) 


] + F ( 


(P.D- 


(32) 


The equations for calculating the points at the right boundary, in cross section X = Xd, are: 

- the downstream boundary condition: 

Qd = Qd it) — outflow hydrograph 

- the positive characteristic equation: 

Q(p+i,d)~Q(p,d) r I7 r -i Q(p+i,d)-Q(p,d-i) , d r T7 , r i f z (p+i,D)~ z (p.D) 


At 

+[ v Cp,d) + C, 


[V 


P,D) 


Ap,D) J 


(P,D)J 


Z(p,D)-Z(p,D-1) 


■ [ Sf(p,D) + ^e(p,D) 


Ax 

1 4- T/ 2 P 4 (p,d)-' 4 (p,d-i) 1 _ _ r 

)J + V (p,D)' [ Ax \ Re(p,D)- L (p,D) 


Ax 

— 9‘ A(P'D) 

' A (-n,D)~ A (p,D- 


' + B (p,Dy [ V (p,D) + C(p,D)] 


At 


+ 


(33) 


Ax Jz=c 

The computation by explicit scheme progresses step by s.tep; the solution is advanced to row (p+1) along the whole length 
of the river stretch that is analyzed. 


IV. DISCUSSION 

The exact mathematical models used for commutating the 
unsteady flow in rivers constitute difficult mathematical 
problems because of the nonlinear nature of the two 
simultaneous partial differential equations, as has already 
been mentioned. Since numerical techniques are too time 
consuming, manual computations are practically 
impossible. However, as it is sometimes far too expensive 
to use digital computers, simplified methods called 
“hydraulic-hydrologic” methods are still preferred in certain 
cases. 

The two partial differential equation (St.Venant’s 
equations), are simplified by omitting one of the two 
equations (momentum equation), or by omitting certain 
terms in this equation. The simplified equations can be 


solved In several different ways, including the numerical 
integration by using manual successive approximation 
procedures. Numerous hydraulic - hydrologic methods of 
flood routing have been developed up to now. Descriptions 
of these methods is recommended for further research. By 
using this method, the outflow hydrograph is obtained as a 
result of calculation. The exactness of the quasi-steady 
method depends on the quality of the storage curve. This 
curve can be defined through hydraulic computations 
(computation of the water levels and of the corresponding 
storage of a river channel), or by balancing the differences 
between the inflow into a river reach and the outflow from 
the reach, for the flood waves registered in the past. The 
rate of storage can be plotted versus discharge. 
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LIST OF DIAGRAMS 



Fig. IB: The net of characteristics 
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Fig. 2: Boundary Conditions 
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